% kpst_graphs.m
%
%  Total versions of the KPST graphs (not per capita or scaled). Their Figures 4 and 5
%   Kelly, B., Papanikolaou, D., Seru, A. and Taddy, M., 2021 AERInsights
%
%   Population is from Statista (Gapminder): https://www.statista.com/statistics/1067138/population-united-states-historical/

clear all; close all;
diarychad('kpst_graphs');
definecolors;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE 2 -- Aggregate Breakthrough patents
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=readtable('Kellyetal2021/time_series_measures_compare.csv'); % Starts in 1836

pop=T.Population_Statista;  % In thousands
Breakthrough = T.Breakthrough_pc.*pop/1000; % T.Breakthrough_pc seems to be per million population in CSV

cshow(' ',[T.year Breakthrough],'%8.0f %20.0f','Year BreakthroughPatents')

startyear = 1900;
t0=find(T.year==startyear);

figure(1); figsetup; makefigwide
plot(T.year(t0:end),Breakthrough(t0:end),'-','Color',myblue,'LineWidth',4);
chadfig('','',1,0);
nums=(0:5000:30000)';
relabelaxis(nums,num2str(nums),'y')
%print kpst_breakthrough.eps

% Log scale
figure(2); figsetup; makefigwide
plot(T.year(t0:end),log(Breakthrough(t0:end)),'-','Color',myblue,'LineWidth',4);
chadfig2('','ratio scale',1,0);
nums=[800 1600 3200 6400 12800 25600]';
relabelaxis(log(nums),num2str(nums),'y')
ax=axis; ax(4)=log(1.05*max(Breakthrough)); axis(ax);
print('-depsc','kpst_breakthrough');

% Growth rates
growthyrs=[1900 1980
           1900 2002
           1950 2002];
yr0=T.year(1)-1;
for t=1:size(growthyrs,1);
    gr(:,t)=100*growthrate(Breakthrough,growthyrs(t,:)-yr0);
end;

disp ' '; disp ' ';
disp 'Growth Rates for KPST Aggregate Data'
cshow(' ',gr,'%12.2f','1900-1980 1900-2002 1950-2002')




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIGURE Industry results
%   Use Breakthrough above to get the "Other" residual category
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tind=readtable('Kellyetal2021/Fig_IndInnovationLR.csv'); % Starts in 1840
indcode=Tind.Properties.VariableNames;
indcode=indcode(2:end); 
FirstYear=Tind.year(1);
LastYear=Tind.year(end);
% Align pop and Industry data to 1840
poporig=pop;
t0pop=find(T.year==FirstYear); % From Breakthrough table
tTpop=find(T.year==LastYear);  % From Breakthrough table
KPST_Industry_pc=Tind{:,2:end};  % The industry data in array form
KPST_Industry = mult(KPST_Industry_pc,pop(t0pop:tTpop))/1000; % Aggregate, not per capita
            % Not clear the /1000 is the right scaling factor here, so initial units maybe scaled?

% Get a TotalBreakthrough that aligns with the Industry years also
TotalBreakthrough=Breakthrough(t0pop:tTpop);
TotalIndustry=sum(KPST_Industry')';
KPST_Industry(:,17) = TotalBreakthrough - TotalIndustry;  % Other = residual between Breakthrough and sum(industry)

% Add the industry labels
name{1} = 'Agriculture, Food';
name{2} = 'Mining, Extraction';
name{3} = 'Utilities';
name{4} = 'Construction';
name{5} = 'Furniture, Textiles, Apparel';
name{6} = 'Wood, Paper, Printing';
name{7} = 'Petroleum, Coal';
name{8} = 'Chemical Manufacturing';
name{9} = 'Plastics, Rubber';
name{10} = 'Mineral Processing';
name{11} = 'Metal Manufacturing';
name{12} = 'Machinery Manufacturing';
name{13} = 'Computers, Electronics';
name{14} = 'Electrical Equipment';
name{15} = 'Transportation Equipment';
name{16} = 'Medical Equipment';
name{17} = 'Other'
disp ' '; disp ' ';
disp 'KPST (2021 AERInsights) Industry Data on Breakthrough Patents (aggregate, not per capita)'
yrstr='1900 1930 1950 1980 2000';
yrs=[1900 1930 1950 1980 2000];
yr0=FirstYear-1;
cshow(name,KPST_Industry(yrs-yr0,:)','%10.1f',yrstr);
Total=sum(KPST_Industry')'; % Now includes Other so should match Breakthrough
cshow('...Total                 ',Total(yrs-yr0)','%10.1f',[],'nonee')


% Growth rates
growthyrs=[1900 1980
           1900 LastYear
           1950 LastYear];

clear gr;
for t=1:size(growthyrs,1);
    gr(:,t)=100*growthrate(KPST_Industry,growthyrs(t,:)-yr0);
end;

disp ' '; disp ' ';
disp 'Growth Rates for KPST Industry Data'
cshow(name,gr,'%12.2f','1900-1980 1900-2002 1950-2002','latex')


% Check that sum of industry is total Breakthrough?


% Aggregate, Agg-Computers, Computers
data=[Total  Total-KPST_Industry(:,13) KPST_Industry(:,13)];
names_summary=[{'Total'} {'Total (excluding computers)'} {'Computers, Electronics'}];
clear gr;
for t=1:size(growthyrs,1);
    gr(:,t)=100*growthrate(data,growthyrs(t,:)-yr0);
end;

disp ' '; disp ' ';
disp 'Growth Rates excluding computers/electronics'
cshow(names_summary,gr,'%12.2f','1900-1980 1900-2002 1950-2002','latex')


diary off;
